#R version 4.4.3 (2025-02-28 ucrt)
#Platform: x86_64-w64-mingw32/x64 (64-bit)
#Running under: Windows 11 x64 

library(dplyr) #1.1.4
library(ggplot2) #3.5.1
library(readr) #2.1.5
library(haven) #2.5.4
library(estimatr) #1.0.6
library(broom.mixed) #0.2.9.6
library(lme4) #1.1-37
library(countrycode) #1.6.1
library(forcats) #1.0.0
library(stargazer) #5.2.3
library(tibble) #3.2.1
library(modelsummary) #2.3.0

`%notin%` <- function(x,y) !(x %in% y) 

setwd("C:/Users/selfd/Documents/autopartyinstitutionalization/legacy/cps/replication")

#load data and set up variables for analysis####
gwf <- read_csv("data/gwf_panel.csv") |>  
  select(countryname, year, v2xps_party, v2x_polyarchy, gwf_duration, gwf_regime, gwf_prior, prior_pi, pi1, pi2, pi3, pi4, asp_dummy, ruling, lgdp, gwf_pduration, gwf_number, fsu, ussr_sat, regime_elecs, av_dm) |>  
  mutate(ccode = countrycode(countryname, origin = "country.name", destination = "cown")) |> 
  mutate_at(vars(pi1, pi2, pi3, pi4), ~ifelse(gwf_prior %in% c("Personal", "Military", "Party", "Monarchy") & is.nan(.) | gwf_prior %in% c("Personal", "Military", "Party", "Monarchy") & is.na(.), 0, .)) |>  
  filter(gwf_prior %in% c("Party", "Military", "Personal") & gwf_regime == "democracy") |>  
  mutate(asp_dummy = ifelse(is.na(asp_dummy), 0, asp_dummy)) |>  
  mutate(gwf_prior = factor(gwf_prior, levels = c("Party", "Personal", "Military"))) |>  
  mutate(elections = ifelse(regime_elecs >= 3, "Third-plus", regime_elecs),  elections = recode(elections, `0` = "None", `1` = "First", `2` = "Second"), 
         elections = factor(elections, levels = c("None", "First", "Second", "Third-plus"))) 

mp <- lm_robust(v2xps_party ~ gwf_prior + pi1 + prior_pi + asp_dummy + v2x_polyarchy + lgdp + gwf_pduration + gwf_number + fsu + ussr_sat + regime_elecs + av_dm, data = gwf, se_type = "HC1") |>
  tidy() |>
  mutate(Sample = "Pooled") |>
  filter(term %in% c("pi1", "prior_pi", "asp_dummy", "gwf_priorParty", "gwf_priorPersonal", "gwf_priorMilitary"))

df1 <- filter(gwf, elections == "None")

m0 <- lm_robust(v2xps_party ~ gwf_prior + pi1 + prior_pi + asp_dummy + v2x_polyarchy + lgdp + gwf_pduration + gwf_number + fsu + ussr_sat + av_dm, data = df1, se_type = "HC1") |>
  tidy() |>
  mutate(Sample = "None") |>
  filter(term %in% c("pi1", "prior_pi", "asp_dummy", "gwf_priorParty", "gwf_priorPersonal", "gwf_priorMilitary"))

df1 <- filter(gwf, elections == "First")

m1 <- lm_robust(v2xps_party ~ gwf_prior + pi1 + prior_pi + asp_dummy + v2x_polyarchy + lgdp + gwf_pduration + gwf_number + fsu + ussr_sat + av_dm, data = df1, se_type = "HC1") |>
  tidy() |>
  mutate(Sample = "First") |>
  filter(term %in% c("pi1", "prior_pi", "asp_dummy", "gwf_priorParty", "gwf_priorPersonal", "gwf_priorMilitary"))

df1 <- filter(gwf, elections == "Second")

m2 <- lm_robust(v2xps_party ~ gwf_prior + pi1 + prior_pi + asp_dummy + v2x_polyarchy + lgdp + gwf_pduration + gwf_number + fsu + ussr_sat + av_dm, data = df1, se_type = "HC1") |>
  tidy() |>
  mutate(Sample = "Second") |>
  filter(term %in% c("pi1", "prior_pi", "asp_dummy", "gwf_priorParty", "gwf_priorPersonal", "gwf_priorMilitary"))

df1 <- filter(gwf, elections == "Third-plus")

m3 <- lm_robust(v2xps_party ~ gwf_prior + pi1 + prior_pi + asp_dummy + v2x_polyarchy + lgdp + gwf_pduration + gwf_number + fsu + ussr_sat + av_dm, data = df1, se_type = "HC1") |>
  tidy() |>
  mutate(Sample = "Third-plus") |>
  filter(term %in% c("pi1", "prior_pi", "asp_dummy", "gwf_priorParty", "gwf_priorPersonal", "gwf_priorMilitary"))

#fig1 left hand side
rbind(m0, m1, m2, m3, mp) |>
  mutate(Sample = factor(Sample, levels = c("None", "First", "Second", "Third-plus", "Pooled"))) |>
  mutate(term = recode(term, "pi1" = "Incumbent PI", "prior_pi" = "Prior PI", "asp_dummy" = "ASP", 
                       "gwf_priorParty" = "Party", "gwf_priorPersonal" = "Personal", "gwf_priorMilitary" = "Military")) |>
  filter(term %in% c("Military", "Personal")) |>
  ggplot(aes(x = Sample, y = estimate, ymin = conf.low, ymax = conf.high, color = term, shape = term)) +
  geom_pointrange(show.legend = TRUE, position = position_jitterdodge(dodge.width = 0.75, jitter.width = 0.2), size = 0.5) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  geom_line(aes(y = estimate - 100), linetype = "solid", color = NA) +  # Setting color to NA to remove legend
  geom_point(aes(y = estimate - 100), shape = NA, color = NA) +  # Setting color and shape to NA to remove legend
  theme_bw() +
  xlab("") +
  ylab("") +
  ylim(-0.55, 0.55) +
  scale_color_grey(name = "Variable", start = 0, end = 0.8) +
  scale_shape_manual("Variable", values = c(15, 16, 17)) + 
  theme(legend.position = "bottom")  +
  ggtitle("OLS Models") +
  theme(plot.title = element_text(hjust = 0.5))

#save to file
#ggsave(filename = "C:/Users/selfd/Documents/autopartyinstitutionalization/legacy/plots/appendix/apb_f1_1.jpg", dpi = 500, width = 5, height = 5)

#RE models#####
gwf8 <- gwf |> select(v2xps_party, gwf_prior, pi1, prior_pi, asp_dummy, v2x_polyarchy, lgdp, gwf_pduration, gwf_number, fsu, ussr_sat, regime_elecs, av_dm, elections, ccode) |> na.omit()

m1 <- lmer(v2xps_party ~ gwf_prior*elections + pi1 + prior_pi + asp_dummy + gwf_pduration + v2x_polyarchy + lgdp + gwf_number + fsu + ussr_sat + regime_elecs + av_dm + (1 | ccode), data = gwf8, REML=FALSE)
m2 <- lmer(v2xps_party ~ gwf_prior*relevel(elections, ref = "First") + pi1 + prior_pi + asp_dummy + gwf_pduration + v2x_polyarchy + lgdp + gwf_number + fsu + ussr_sat + regime_elecs + av_dm + (1 | ccode), data = gwf8, REML=FALSE)
m3 <- lmer(v2xps_party ~ gwf_prior*relevel(elections, ref = "Second") + pi1 + prior_pi + asp_dummy + gwf_pduration + v2x_polyarchy + lgdp + gwf_number + fsu + ussr_sat + regime_elecs + av_dm + (1 | ccode), data = gwf8, REML=FALSE)
m4 <- lmer(v2xps_party ~ gwf_prior*relevel(elections, ref = "Third-plus") + pi1 + prior_pi + asp_dummy + gwf_pduration + v2x_polyarchy + lgdp + gwf_number + fsu + ussr_sat + regime_elecs + av_dm + (1 | ccode), data = gwf8, REML=FALSE)

#get pi1 effect for each round of elections
tm1 <- tidy(m1, conf.int = TRUE) |> filter(term %in% c("gwf_priorPersonal", "gwf_priorMilitary")) |> mutate(Sample = "None")

tm2 <- tidy(m2, conf.int = TRUE) |> filter(term %in% c("gwf_priorPersonal", "gwf_priorMilitary")) |> mutate(Sample = "First")

tm3 <- tidy(m3, conf.int = TRUE) |> filter(term %in% c("gwf_priorPersonal", "gwf_priorMilitary")) |> mutate(Sample = "Second")

tm4 <- tidy(m4, conf.int = TRUE) |> filter(term %in% c("gwf_priorPersonal", "gwf_priorMilitary")) |> mutate(Sample = "Third-plus")

rbind(tm1, tm2, tm3, tm4) |> 
  mutate(term = case_when(term == "gwf_priorPersonal" ~ "Personal", term == "gwf_priorMilitary" ~ "Military"),
         Sample = fct_relevel(Sample, c("None", "First", "Second", "Third-plus"))) |> 
  ggplot(aes(x = Sample, y = estimate, ymin = conf.low, ymax = conf.high, color = term, shape = term)) +
  geom_pointrange(show.legend = T, position = position_jitterdodge(dodge.width = 0.75, jitter.width = 0.2), size = 0.5) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  geom_line(aes(y = estimate - 100), linetype = "solid", color = NA) +  # Setting color to NA to remove legend
  geom_point(aes(y = estimate - 100), shape = NA, color = NA) +  # Setting color and shape to NA to remove legend
  theme_bw() +
  xlab("") +
  ylab("") +
  ylim(-0.55, 0.55) +
  scale_color_grey(name = "Variable", start = 0, end = 0.8) +
  scale_shape_manual("Variable", values = c(15, 16, 17)) + 
  theme(legend.position = "bottom")  +
  ggtitle("Random Effects Models") +
  theme(plot.title = element_text(hjust = 0.5))

#save to file
#ggsave(filename = "C:/Users/selfd/Documents/autopartyinstitutionalization/legacy/plots/appendix/apb_f1_2.jpg", dpi = 500, width = 5, height = 5)

#Table 1, OLS Models------
#Ran in stata to extract PCSE using an unbalanced panel
#See model1.1, model1.2, model1.3 model1.4 and model1.5 with output from regression_results.tex

#Table 2, RE Models ----
gwf <- read_csv("data/gwf_panel.csv") %>% 
  select(countryname, year, v2xps_party, v2x_polyarchy, gwf_duration, gwf_regime, gwf_prior, prior_pi, pi1, pi2, pi3, pi4, asp_dummy, ruling, lgdp, gwf_pduration, gwf_number, fsu, ussr_sat, regime_elecs, av_dm) %>% 
  mutate(ccode = countrycode(countryname, origin = "country.name", destination = "cown")) |> 
  mutate_at(vars(pi1, pi2, pi3, pi4), ~ifelse(gwf_prior %in% c("Personal", "Military", "Party", "Monarchy") & is.nan(.) | gwf_prior %in% c("Personal", "Military", "Party", "Monarchy") & is.na(.), 0, .)) %>% 
  filter(gwf_prior %in% c("Party", "Military", "Personal") & gwf_regime == "democracy") %>% 
  mutate(asp_dummy = ifelse(is.na(asp_dummy), 0, asp_dummy)) %>% 
  mutate(gwf_prior = factor(gwf_prior, levels = c("Party", "Personal", "Military"))) %>% 
  mutate(elections = ifelse(regime_elecs >= 3, "Third-plus", regime_elecs),  elections = recode(elections, `0` = "None", `1` = "First", `2` = "Second"), 
         elections = factor(elections, levels = c("None", "First", "Second", "Third-plus"))) 

gwf8 <- gwf |> select(v2xps_party, gwf_prior, pi1, prior_pi, asp_dummy, v2x_polyarchy, lgdp, gwf_pduration, gwf_number, fsu, ussr_sat, regime_elecs, av_dm, elections, ccode) |> na.omit()

m1 <- lmer(v2xps_party ~ gwf_prior + pi1*elections + prior_pi + asp_dummy + (1 | ccode), data = gwf8, REML=FALSE)
m2 <- lmer(v2xps_party ~ gwf_prior + pi1*relevel(elections, ref = "First") + prior_pi + asp_dummy + (1 | ccode), data = gwf8, REML=FALSE)
m3 <- lmer(v2xps_party ~ gwf_prior + pi1*relevel(elections, ref = "Second") + prior_pi + asp_dummy + (1 | ccode), data = gwf8, REML=FALSE)
m4 <- lmer(v2xps_party ~ gwf_prior + pi1*relevel(elections, ref = "Third-plus") + prior_pi + asp_dummy + + av_dm + (1 | ccode), data = gwf8, REML=FALSE)

#get pi1 effect for each round of elections
tm1 <- tidy(m1, conf.int = TRUE) |> filter(term == "pi1") |> mutate(Sample = "None")

tm2 <- tidy(m2, conf.int = TRUE) |> filter(term == "pi1") |> mutate(Sample = "First")

tm3 <- tidy(m3, conf.int = TRUE) |> filter(term == "pi1") |> mutate(Sample = "Second")

tm4 <- tidy(m4, conf.int = TRUE) |> filter(term == "pi1") |> mutate(Sample = "Third-plus")

#now do it for prior_pi
m5 <- lmer(v2xps_party ~ gwf_prior + prior_pi*elections + pi1 + asp_dummy + (1 | ccode), data = gwf8, REML=FALSE)
m6 <- lmer(v2xps_party ~ gwf_prior + prior_pi*relevel(elections, ref = "First") + pi1 + asp_dummy + (1 | ccode), data = gwf8, REML=FALSE)
m7 <- lmer(v2xps_party ~ gwf_prior + prior_pi*relevel(elections, ref = "Second") + pi1 + asp_dummy + (1 | ccode), data = gwf8, REML=FALSE)
m8 <- lmer(v2xps_party ~ gwf_prior + prior_pi*relevel(elections, ref = "Third-plus") + pi1 + asp_dummy + (1 | ccode), data = gwf8, REML=FALSE)

#get pi1 effect for each round of elections
tm5 <- tidy(m5, conf.int = TRUE) |> filter(term == "prior_pi") |> mutate(Sample = "None")

tm6 <- tidy(m6, conf.int = TRUE) |> filter(term == "prior_pi") |> mutate(Sample = "First")

tm7 <- tidy(m7, conf.int = TRUE) |> filter(term == "prior_pi") |> mutate(Sample = "Second")

tm8 <- tidy(m8, conf.int = TRUE) |> filter(term == "prior_pi") |> mutate(Sample = "Third-plus")

#lastly we do it for asp
m9 <- lmer(v2xps_party ~ gwf_prior + asp_dummy*elections + pi1 + prior_pi + (1 | ccode), data = gwf8, REML=FALSE)
m10 <- lmer(v2xps_party ~ gwf_prior + asp_dummy*relevel(elections, ref = "First") + pi1 + prior_pi + (1 | ccode), data = gwf8, REML=FALSE)
m11 <- lmer(v2xps_party ~ gwf_prior + asp_dummy*relevel(elections, ref = "Second") + pi1 + prior_pi + (1 | ccode), data = gwf8, REML=FALSE)
m12 <- lmer(v2xps_party ~ gwf_prior + asp_dummy*relevel(elections, ref = "Third-plus") + pi1 + prior_pi + (1 | ccode), data = gwf8, REML=FALSE)

#get pi1 effect for each round of elections
tm9 <- tidy(m9, conf.int = TRUE) |> filter(term == "asp_dummy") |> mutate(Sample = "None")

tm10 <- tidy(m10, conf.int = TRUE) |> filter(term == "asp_dummy") |> mutate(Sample = "First")

tm11 <- tidy(m11, conf.int = TRUE) |> filter(term == "asp_dummy") |> mutate(Sample = "Second")

tm12 <- tidy(m12, conf.int = TRUE) |> filter(term == "asp_dummy") |> mutate(Sample = "Third-plus")

modellist <- list("Pooled - PI" = m1, 
                  "First - PI" = m2, 
                  "Second - PI" = m3, 
                  "Third+ - PI" = m4, 
                  "Pooled - Prior" = m5, 
                  "First - Prior" = m6, 
                  "Second - Prior" = m7, 
                  "Third+ - Prior" = m8, 
                  "Pooled - ASP" = m9, 
                  "First - ASP" = m10, 
                  "Second - ASP" = m11, 
                  "Third+ - ASP" = m12) 

modelsummary(modellist, stars = TRUE, gof_omit = "R2|AIC|RM|F|Log.Lik", coef_map = c("pi1" = "Incumbent PI", "prior_pi" = "Prior PI", "asp_dummy" = "ASP",
                                                                                     "reactive" = "Reactive ASP",
                                                                                     "gwf_priorPersonal" = "Personalist", "gwf_priorMilitary" = "Military",
                                                                                     "gwf_pduration" = "Prior Duration", "regime_elecs" = "Elections", 
                                                                                     "SD (Intercept ccode)" = "SD: Country", "SD (Observations)" = "SD: Residuals"), output = "latex")



#Figure 3, ols and re model ----
gwf <- read_csv("data/gwf_panel.csv") |>  
  select(countryname, year, v2xps_party, v2x_polyarchy, gwf_duration, gwf_regime, gwf_prior, prior_pi, pi1, pi2, pi3, pi4, asp_dummy, ruling, lgdp, gwf_pduration, gwf_number, fsu, ussr_sat, regime_elecs, av_dm) |>  
  mutate(ccode = countrycode(countryname, origin = "country.name", destination = "cown"), reactive = case_when(ruling == 1 ~ 0, ruling == 0 ~ 1)) |> 
  mutate_at(vars(pi1, pi2, pi3, pi4), ~ifelse(gwf_prior %in% c("Personal", "Military", "Party", "Monarchy") & is.nan(.) | gwf_prior %in% c("Personal", "Military", "Party", "Monarchy") & is.na(.), 0, .)) |>  
  filter(gwf_prior %in% c("Party", "Military", "Personal") & gwf_regime == "democracy") |>  
  mutate(asp_dummy = ifelse(is.na(asp_dummy), 0, asp_dummy)) |>  
  mutate(gwf_prior = factor(gwf_prior, levels = c("Party", "Personal", "Military"))) |>  
  mutate(elections = ifelse(regime_elecs >= 3, "Third-plus", regime_elecs),  elections = recode(elections, `0` = "None", `1` = "First", `2` = "Second"), 
         elections = factor(elections, levels = c("None", "First", "Second", "Third-plus"))) 

mp <- lm_robust(v2xps_party ~ gwf_prior + pi1 + prior_pi + reactive + v2x_polyarchy + lgdp + gwf_pduration + gwf_number + fsu + ussr_sat + regime_elecs + av_dm, data = gwf, se_type = "HC1")
df1 <- filter(gwf, elections == "None")

m0 <- lm_robust(v2xps_party ~ gwf_prior + pi1 + prior_pi + reactive + v2x_polyarchy + lgdp + gwf_pduration + gwf_number + fsu + ussr_sat + av_dm, data = df1, se_type = "HC1") 
df1 <- filter(gwf, elections == "First")

m1 <- lm_robust(v2xps_party ~ gwf_prior + pi1 + prior_pi + reactive + v2x_polyarchy + lgdp + gwf_pduration + gwf_number + fsu + ussr_sat + av_dm, data = df1, se_type = "HC1") 
df1 <- filter(gwf, elections == "Second")

m2 <- lm_robust(v2xps_party ~ gwf_prior + pi1 + prior_pi + reactive + v2x_polyarchy + lgdp + gwf_pduration + gwf_number + fsu + ussr_sat + av_dm, data = df1, se_type = "HC1") 
df1 <- filter(gwf, elections == "Third-plus")

m3 <- lm_robust(v2xps_party ~ gwf_prior + pi1 + prior_pi + reactive + v2x_polyarchy + lgdp + gwf_pduration + gwf_number + fsu + ussr_sat + av_dm, data = df1, se_type = "HC1") 

gwf8 <- gwf |> select(v2xps_party, gwf_prior, pi1, prior_pi, reactive, v2x_polyarchy, lgdp, gwf_pduration, gwf_number, fsu, ussr_sat, regime_elecs, av_dm, elections, ccode) |> na.omit()

m4 <- lmer(v2xps_party ~ gwf_prior + pi1*elections + prior_pi + reactive + gwf_pduration + v2x_polyarchy + lgdp + gwf_number + fsu + ussr_sat + regime_elecs + av_dm + (1 | ccode), data = gwf8, REML=FALSE)
m5 <- lmer(v2xps_party ~ gwf_prior + prior_pi*elections + pi1 + reactive + gwf_pduration + v2x_polyarchy + lgdp + gwf_number + fsu + ussr_sat + regime_elecs + av_dm + (1 | ccode), data = gwf8, REML=FALSE)
m6 <- lmer(v2xps_party ~ gwf_prior + reactive*elections + prior_pi + pi1 + gwf_pduration + v2x_polyarchy + lgdp + gwf_number + fsu + ussr_sat + regime_elecs + av_dm + (1 | ccode), data = gwf8, REML=FALSE)


modellist <- list("OLS - None" <- m0,
                  "OLS - First" <- m1,
                  "OLS - Second" <- m2,
                  "OLS - Third + " <- m3,
                  "OLS - Pooled" <- mp,
                  "ML - Incumbent PI" <- m4,
                  "ML - Prior PI" <- m5,
                  "ML - Reactive ASP" <- m6)

modelsummary(modellist, stars = TRUE, gof_omit = "R2|AIC|RM|F|Log.Lik", output = "latex")


#Interaction models Fiugres 3, 4, and 5-----
setwd("legacy")
df1 <- read_csv("data/gwf_panel.csv") |>  
  dplyr::select(countryname, year, v2xps_party, v2x_polyarchy, gwf_duration, gwf_regime, gwf_prior, prior_pi, pi1, pi2, pi3, pi4, asp_dummy, ruling, lgdp, gwf_pduration, gwf_number, fsu, ussr_sat, regime_elecs, av_dm) |>  
  mutate(ccode = countrycode(countryname, origin = "country.name", destination = "cown"), 
         reactive = case_when(ruling == 1 ~ 0, ruling == 0 ~ 1)) |> 
  mutate_at(vars(pi1, pi2, pi3, pi4), ~ifelse(gwf_prior %in% c("Personal", "Military", "Party", "Monarchy") & is.nan(.) | gwf_prior %in% c("Personal", "Military", "Party", "Monarchy") & is.na(.), 0, .)) |>  
  filter(gwf_prior %in% c("Party", "Military", "Personal") & gwf_regime == "democracy") |>  
  mutate(asp_dummy = ifelse(is.na(asp_dummy), 0, asp_dummy)) |>  
  mutate(gwf_prior = factor(gwf_prior, levels = c("Party", "Personal", "Military"))) |>  
  mutate(elections = ifelse(regime_elecs >= 3, "Third-plus", regime_elecs),  elections = recode(elections, `0` = "None", `1` = "First", `2` = "Second"), 
         elections = factor(elections, levels = c("None", "First", "Second", "Third-plus"))) |> 
  select(v2xps_party, gwf_prior, reactive, pi1, pi4, prior_pi, v2x_polyarchy, lgdp, fsu, ussr_sat, av_dm) |> 
  na.omit()
#Figure 3
im1 <- lm(v2xps_party ~ gwf_prior + reactive*pi1 + prior_pi + v2x_polyarchy + lgdp + fsu + ussr_sat + av_dm, data = df1)
im2 <- lm(v2xps_party ~ gwf_prior + reactive*prior_pi + pi1 + v2x_polyarchy + lgdp + fsu + ussr_sat + av_dm, data = df1)
im3 <- lm(v2xps_party ~ gwf_prior + reactive*pi4 + prior_pi + v2x_polyarchy + lgdp + fsu + ussr_sat + av_dm, data = df1)

modellist <- list("Figure 3" <- im1,
                  "Figure 4" <- im2,
                  "Figure 5" <- im3) 
modelsummary(modellist, stars = TRUE, gof_omit = "R2|AIC|RM|F|Log.Lik", output = "latex")



